- Vincent Fugère (il/lui)
- prof. en sciences de l’environnement à l’UQTR
- recherche en écologie des eaux douces au RIVE
- science des données & environnement
02-2022
dplyr, ggplot2lm()Nous allons réutiliser les ensembles de données que vous avez utilisé lors de deux atelier précédents: gapminder, palmer pinguins, starwars. Nous devons également charger les librairies dplyr et ggplot2.
library(dplyr)
library(ggplot2)
library(palmerpenguins)
library(gapminder)
data("starwars")
data("penguins")
data("gapminder")
# #au besoin
# install.packages(c("penguins","gapminder"))
head(penguins)
## # A tibble: 6 × 8 ## species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex ## <fct> <fct> <dbl> <dbl> <int> <int> <fct> ## 1 Adelie Torge… 39.1 18.7 181 3750 male ## 2 Adelie Torge… 39.5 17.4 186 3800 fema… ## 3 Adelie Torge… 40.3 18 195 3250 fema… ## 4 Adelie Torge… NA NA NA NA <NA> ## 5 Adelie Torge… 36.7 19.3 193 3450 fema… ## 6 Adelie Torge… 39.3 20.6 190 3650 male ## # … with 1 more variable: year <int>
Identifiez les variables qualitatives et quantitatives de l’ensemble de données starwars
Identifiez les variables qualitatives et quantitatives de l’ensemble de données starwars
head(starwars)
glimpse(starwars)
## Rows: 87 ## Columns: 14 ## $ name <chr> "Luke Skywalker", "C-3PO", "R2-D2", "Darth Vader", "Leia Or… ## $ height <int> 172, 167, 96, 202, 150, 178, 165, 97, 183, 182, 188, 180, 2… ## $ mass <dbl> 77.0, 75.0, 32.0, 136.0, 49.0, 120.0, 75.0, 32.0, 84.0, 77.… ## $ hair_color <chr> "blond", NA, NA, "none", "brown", "brown, grey", "brown", N… ## $ skin_color <chr> "fair", "gold", "white, blue", "white", "light", "light", "… ## $ eye_color <chr> "blue", "yellow", "red", "yellow", "brown", "blue", "blue",… ## $ birth_year <dbl> 19.0, 112.0, 33.0, 41.9, 19.0, 52.0, 47.0, NA, 24.0, 57.0, … ## $ sex <chr> "male", "none", "none", "male", "female", "male", "female",… ## $ gender <chr> "masculine", "masculine", "masculine", "masculine", "femini… ## $ homeworld <chr> "Tatooine", "Tatooine", "Naboo", "Tatooine", "Alderaan", "T… ## $ species <chr> "Human", "Droid", "Droid", "Human", "Human", "Human", "Huma… ## $ films <list> <"The Empire Strikes Back", "Revenge of the Sith", "Return… ## $ vehicles <list> <"Snowspeeder", "Imperial Speeder Bike">, <>, <>, <>, "Imp… ## $ starships <list> <"X-wing", "Imperial shuttle">, <>, <>, "TIE Advanced x1",…
\[ y = mx + b \] \[ y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_px_{ip} + \epsilon_i \] \[ \epsilon \sim N(0,\sigma^2) \]
hist(rnorm(10000))
penguins %>% filter(species == 'Adelie') %>% group_by(sex) %>% summarize(masse.moy = mean(body_mass_g))
## # A tibble: 3 × 2 ## sex masse.moy ## <fct> <dbl> ## 1 female 3369. ## 2 male 4043. ## 3 <NA> NA
penguins %>% filter(species == 'Adelie', sex != 'NA') %>% ggplot(aes(x=body_mass_g, colour=sex)) + geom_density()
Il semble y avoir une différence de masse entre les sexes. La fonction t.test nous permet de valider cette hypothèse.
adelie <- penguins %>% filter(species == 'Adelie', sex != 'NA') t.test(body_mass_g ~ sex, data = adelie)
## ## Welch Two Sample t-test ## ## data: body_mass_g by sex ## t = -13.126, df = 135.69, p-value < 2.2e-16 ## alternative hypothesis: true difference in means between group female and group male is not equal to 0 ## 95 percent confidence interval: ## -776.3012 -573.0139 ## sample estimates: ## mean in group female mean in group male ## 3368.836 4043.493
Il semble y avoir une différence de masse entre les sexes. La fonction t.test nous permet de valider cette hypothèse.
adelie <- penguins %>% filter(species == 'Adelie', sex != 'NA') t.test(body_mass_g ~ sex, data = adelie)
## ## Welch Two Sample t-test ## ## data: body_mass_g by sex ## t = -13.126, df = 135.69, p-value < 2.2e-16 ## alternative hypothesis: true difference in means between group female and group male is not equal to 0 ## 95 percent confidence interval: ## -776.3012 -573.0139 ## sample estimates: ## mean in group female mean in group male ## 3368.836 4043.493
format(2.2e-16, scientific = F)
## [1] "0.00000000000000022"
La fonction lm fonctionne pour tous les types de modèles linéaires (incluant le test de t).
modele <- lm(body_mass_g ~ sex, data = adelie) summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ sex, data = adelie) ## ## Residuals: ## Min 1Q Median 3Q Max ## -718.49 -218.84 -18.84 225.00 731.51 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3368.84 36.34 92.69 <2e-16 *** ## sexmale 674.66 51.40 13.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 310.5 on 144 degrees of freedom ## Multiple R-squared: 0.5447, Adjusted R-squared: 0.5416 ## F-statistic: 172.3 on 1 and 144 DF, p-value: < 2.2e-16
\[ Y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \] Quelles valeurs la variable “x” (sex) peut-elle prendre?
summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ sex, data = adelie) ## ## Residuals: ## Min 1Q Median 3Q Max ## -718.49 -218.84 -18.84 225.00 731.51 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3368.84 36.34 92.69 <2e-16 *** ## sexmale 674.66 51.40 13.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 310.5 on 144 degrees of freedom ## Multiple R-squared: 0.5447, Adjusted R-squared: 0.5416 ## F-statistic: 172.3 on 1 and 144 DF, p-value: < 2.2e-16
Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo? On peut extraire les données pertinentes comme suit:
nab_tat <- starwars %>% filter(homeworld %in% c('Naboo','Tatooine'))
ggplot(aes(x=homeworld,y=height),data=nab_tat) + geom_boxplot()
Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo?
t.test(height ~ homeworld, data = nab_tat)
## ## Welch Two Sample t-test ## ## data: height by homeworld ## t = 0.42386, df = 18.947, p-value = 0.6764 ## alternative hypothesis: true difference in means between group Naboo and group Tatooine is not equal to 0 ## 95 percent confidence interval: ## -22.27276 33.58185 ## sample estimates: ## mean in group Naboo mean in group Tatooine ## 175.4545 169.8000
Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo?
modele <- lm(height ~ homeworld, data = nab_tat) summary(modele)
## ## Call: ## lm(formula = height ~ homeworld, data = nab_tat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -79.455 -6.800 7.545 13.200 48.545 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 175.455 9.276 18.914 8.76e-14 *** ## homeworldTatooine -5.655 13.443 -0.421 0.679 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 30.77 on 19 degrees of freedom ## Multiple R-squared: 0.009227, Adjusted R-squared: -0.04292 ## F-statistic: 0.1769 on 1 and 19 DF, p-value: 0.6787
Lorsque notre variable qualitative indépendante comprend plus de deux niveaux, nous pouvons utiliser une analyse de variance (ANOVA) plutôt qu’un test de t. On ajuste le modèle de la même façon avec la fonction lm:
modele <- lm(body_mass_g ~ species, data = penguins) summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ species, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1126.02 -333.09 -33.09 316.91 1223.98 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3700.66 37.62 98.37 <2e-16 *** ## speciesChinstrap 32.43 67.51 0.48 0.631 ## speciesGentoo 1375.35 56.15 24.50 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 462.3 on 339 degrees of freedom ## (2 observations deleted due to missingness) ## Multiple R-squared: 0.6697, Adjusted R-squared: 0.6677 ## F-statistic: 343.6 on 2 and 339 DF, p-value: < 2.2e-16
Les coefficients (béta) du modèle indiquent la différence de masse entre chaque niveau et le niveau de référence. Par défaut, le niveau de référence d’un facteur est celui qui arrive en première position lorsqu’on classe les niveaux en ordre alphabétique. Pour nos manchots, il y a trois espèces: Adelie, Chinstrap, Gentoo. Adelie est donc le niveau de référence.
levels(penguins$species)
## [1] "Adelie" "Chinstrap" "Gentoo"
coef(modele)
## (Intercept) speciesChinstrap speciesGentoo ## 3700.66225 32.42598 1375.35401
masse moyenne Adelie = 3700, Chinstrap = 3700+32, Gentoo = 3700+1375
masse moyenne Adelie = 3700, Chinstrap = 3700+32, Gentoo = 3700+1375
ggplot(aes(x=species,y=body_mass_g),data=penguins) + geom_boxplot()
Souvent, pour une ANOVA, on cherche plutôt à tester si la variable indépendante influence au moins un des niveaux, peu importe lequel, et sans nécessairement comparer les niveaux à un niveau de référence. On veut un tableau d’ANOVA qui ressemble à cela:
Il suffit d’utiliser la fonction anova pour obtenir le tableau d’ANOVA de notre modèle:
anova(modele)
## Analysis of Variance Table ## ## Response: body_mass_g ## Df Sum Sq Mean Sq F value Pr(>F) ## species 2 146864214 73432107 343.63 < 2.2e-16 *** ## Residuals 339 72443483 213698 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
On complémente souvent le test de F de l’ANOVA avec des comparaisons par paires, pour déterminer quels groupes diffèrent les uns des autres. Dans ce cas, y a-t-il une différence entre Chinstrap et Adelie?
ggplot(aes(x=species,y=body_mass_g),data=penguins) + geom_boxplot()
Comparaisons par paires à l’aide de tests de Tukey (“Tukey’s Honest Significant Difference”):
TukeyHSD(aov(modele))
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = modele) ## ## $species ## diff lwr upr p adj ## Chinstrap-Adelie 32.42598 -126.5002 191.3522 0.8806666 ## Gentoo-Adelie 1375.35401 1243.1786 1507.5294 0.0000000 ## Gentoo-Chinstrap 1342.92802 1178.4810 1507.3750 0.0000000
À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?
donnees_2007 <- gapminder %>% filter(year == 2007) ggplot(aes(x=continent,y=lifeExp),data=donnees_2007) + geom_boxplot()
À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?
modele <- lm(lifeExp ~ continent, donnees_2007) anova(modele)
## Analysis of Variance Table ## ## Response: lifeExp ## Df Sum Sq Mean Sq F value Pr(>F) ## continent 4 13060.7 3265.2 59.714 < 2.2e-16 *** ## Residuals 137 7491.2 54.7 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?
TukeyHSD(aov(modele))
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = modele) ## ## $continent ## diff lwr upr p adj ## Americas-Africa 18.802082 13.827042 23.777121 0.0000000 ## Asia-Africa 15.922446 11.372842 20.472051 0.0000000 ## Europe-Africa 22.842562 18.155858 27.529265 0.0000000 ## Oceania-Africa 25.913462 11.183451 40.643472 0.0000305 ## Asia-Americas -2.879635 -8.299767 2.540497 0.5844901 ## Europe-Americas 4.040480 -1.495233 9.576193 0.2630707 ## Oceania-Americas 7.111380 -7.910342 22.133102 0.6862848 ## Europe-Asia 6.920115 1.763372 12.076858 0.0027416 ## Oceania-Asia 9.991015 -4.895221 24.877251 0.3464970 ## Oceania-Europe 3.070900 -11.857807 17.999607 0.9793776
Merveilleux, mais… en regardant ces boîtes à moustaches (boxplot), trouvez-vous qu’il y a quelque chose de louche?
ggplot(aes(x=continent,y=lifeExp),data=donnees_2007) + geom_boxplot()
\[ Y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \]
\[ \epsilon \sim N(0,\sigma^2) \]
Quelques suppositions des modèles linéaires:
Nous pouvons vérifier ces suppositions à l’aide de tests. IMPORTANT: pour ces deux tests, nous voulons que p > 0.05!
#1. normalité des résidus mod.resid <- resid(modele) # Sortir les résidus shapiro.test(mod.resid) # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.96708, p-value = 0.001696
Nous pouvons vérifier ces suppositions à l’aide de tests. IMPORTANT: pour ces deux tests, nous voulons que p > 0.05!
#1. normalité des résidus mod.resid <- resid(modele) # Sortir les résidus shapiro.test(mod.resid) # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.96708, p-value = 0.001696
# 2. homoscédasticité bartlett.test(lifeExp ~ continent, donnees_2007) # Test de Bartlett: p < 0.05 = MAUVAIS!
## ## Bartlett test of homogeneity of variances ## ## data: lifeExp by continent ## Bartlett's K-squared = 45.85, df = 4, p-value = 2.646e-09
Il faut vérifier les suppositions d’un modèle avant de l’interpréter! S’il y a un problème, on ne regarde même pas les résultats.
Solutions potentielles:
Une ANOVA peut inclure plus d’une variables qualitatives, ce qui permet notamment de quantifier les intéractions entre deux ou plusieurs variables. Par exemple, est-ce que l’effet du sexe sur la masse des manchots dépend de l’espèce?
penguins <- tidyr::drop_na(penguins) ggplot(aes(x=species,y=body_mass_g,fill=sex),data=penguins) + geom_boxplot()
Une ANOVA peut inclure plus d’une variables qualitatives, ce qui permet notamment de quantifier les intéractions entre deux ou plusieurs variables. Par exemple, est-ce que l’effet du sexe sur la masse des manchots dépend de l’espèce?
modele <- lm(body_mass_g ~ species + sex, data = penguins) anova(modele)
## Analysis of Variance Table ## ## Response: body_mass_g ## Df Sum Sq Mean Sq F value Pr(>F) ## species 2 145190219 72595110 724.21 < 2.2e-16 *** ## sex 1 37090262 37090262 370.01 < 2.2e-16 *** ## Residuals 329 32979185 100241 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Que manque-t-il?
modele <- lm(body_mass_g ~ species * sex, data = penguins) anova(modele)
## Analysis of Variance Table ## ## Response: body_mass_g ## Df Sum Sq Mean Sq F value Pr(>F) ## species 2 145190219 72595110 758.358 < 2.2e-16 *** ## sex 1 37090262 37090262 387.460 < 2.2e-16 *** ## species:sex 2 1676557 838278 8.757 0.0001973 *** ## Residuals 327 31302628 95727 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
N’oublions pas qu’il s’agit encore du modèle linéaire général:
\[ Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_px_{ip} + \epsilon_i \]
summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ species * sex, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -827.21 -213.97 11.03 206.51 861.03 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3368.84 36.21 93.030 < 2e-16 *** ## speciesChinstrap 158.37 64.24 2.465 0.01420 * ## speciesGentoo 1310.91 54.42 24.088 < 2e-16 *** ## sexmale 674.66 51.21 13.174 < 2e-16 *** ## speciesChinstrap:sexmale -262.89 90.85 -2.894 0.00406 ** ## speciesGentoo:sexmale 130.44 76.44 1.706 0.08886 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 309.4 on 327 degrees of freedom ## Multiple R-squared: 0.8546, Adjusted R-squared: 0.8524 ## F-statistic: 384.3 on 5 and 327 DF, p-value: < 2.2e-16
Chez les personnages de star wars, y a-t-il une relation entre la masse corporelle et la grandeur?
ggplot(aes(x=height,y=mass),data=starwars) + geom_point()
Chez les personnages de star wars, y a-t-il une relation entre la masse corporelle et la grandeur?
ggplot(aes(x=height,y=mass),data=starwars) + geom_point()
Chez les personnages de star wars (autres que Jabba the Hut), y a-t-il une relation entre la masse corporelle et la grandeur?
sw <- starwars %>% filter(mass < 1000) ggplot(aes(x=height,y=mass),data=sw) + geom_point()
Chez les personnages de star wars (autres que Jabba the Hut), y a-t-il une relation entre la masse corporelle et la grandeur?
modele <- lm(mass ~ height, data = sw) summary(modele)
## ## Call: ## lm(formula = mass ~ height, data = sw) ## ## Residuals: ## Min 1Q Median 3Q Max ## -39.382 -8.212 0.211 3.846 57.327 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -32.54076 12.56053 -2.591 0.0122 * ## height 0.62136 0.07073 8.785 4.02e-12 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.14 on 56 degrees of freedom ## Multiple R-squared: 0.5795, Adjusted R-squared: 0.572 ## F-statistic: 77.18 on 1 and 56 DF, p-value: 4.018e-12
Validons les suppositions. Premièrement, la normalité des résidus:
mod.resid <- resid(modele) shapiro.test(mod.resid) # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.912, p-value = 0.0004688
Ça ne passe pas…
On ré-essaie en transformant notre variable dépendante.
sw$log_masse <- log10(sw$mass) modele <- lm(log_masse ~ height, data = sw) mod.resid <- resid(modele) shapiro.test(mod.resid) #yé!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.96702, p-value = 0.1156
Pour l’homogénéité de la variance, il suffit de regarder la distribution des résidus en fonction des valeurs de y estimées par le modèle.
Graphique tiré de ce site web
Pour l’homogénéité de la variance, il suffit de regarder la distribution des résidus en fonction des valeurs de y estimées par le modèle. Ici, on espère observer un nuage de point sans aucune forme triangulaire ou tendance directionnelle évidente.
plot(modele, which = 1)
Tout est bien. On peut interpréter le modèle…
summary(modele)
## ## Call: ## lm(formula = log_masse ~ height, data = sw) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.25043 -0.06977 0.01734 0.05502 0.22060 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.9911872 0.0690142 14.36 <2e-16 *** ## height 0.0048730 0.0003886 12.54 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.1052 on 56 degrees of freedom ## Multiple R-squared: 0.7374, Adjusted R-squared: 0.7327 ## F-statistic: 157.2 on 1 and 56 DF, p-value: < 2.2e-16
Tout est bien. On peut interpréter le modèle et le visualiser.
ggplot(aes(x=height,y=log_masse),data=sw) + geom_point() + geom_smooth(method='lm')
\[ y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \]
Plutôt que d’utiliser ggplot pour ajuster le modèle à nouveau, nous pouvons aussi extraire les coéfficients (béta) de notre modèle et les utiliser pour “dessiner” la droite.
coef(modele)
## (Intercept) height ## 0.991187162 0.004872992
b0 <- coef(modele)[1] b1 <- coef(modele)[2] c(b0,b1)
## (Intercept) height ## 0.991187162 0.004872992
Quelle est l’équation de notre modèle linéaire?
\[ y_i = 0.9912 + 0.0049 * height_i + \epsilon_i \]
ggplot(aes(x=height,y=log_masse),data=sw) + geom_point() + geom_abline(intercept = b0, slope = b1)
À votre tour! Dans les données penguins, y a-t-il une relation entre la longueur du bec (bill_length_mm) et la longueur des nageoires (flipper_length_mm)?
modele <- lm(bill_length_mm ~ flipper_length_mm, penguins) par(mfrow=c(1,2)) plot(modele, which = 1:2)
summary(modele)
## ## Call: ## lm(formula = bill_length_mm ~ flipper_length_mm, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -8.6367 -2.6981 -0.5788 2.0663 19.0953 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -7.21856 3.27175 -2.206 0.028 * ## flipper_length_mm 0.25482 0.01624 15.691 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 4.148 on 331 degrees of freedom ## Multiple R-squared: 0.4265, Adjusted R-squared: 0.4248 ## F-statistic: 246.2 on 1 and 331 DF, p-value: < 2.2e-16
longueur du bec = -7.21856 + 0.25482 * longueur de la nageoire
ggplot(aes(x=flipper_length_mm,y=bill_length_mm),data=penguins) + geom_point() + geom_smooth(method='lm')
Analyse de covariance: 1 variable qualitative, 1 variable quantitative.
ggplot(aes( x=body_mass_g, y=bill_length_mm, colour = species), data = penguins) + geom_point() + geom_smooth(method = "lm", aes(fill = species))
modele.sans.interaction <- lm(bill_length_mm ~ body_mass_g + species, data = penguins) summary(modele.sans.interaction) #equation pour le manchot Gentoo?
## ## Call: ## lm(formula = bill_length_mm ~ body_mass_g + species, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.8291 -1.6728 0.1244 1.5318 9.2904 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 24.908763 1.089730 22.858 < 2e-16 *** ## body_mass_g 0.003755 0.000289 12.991 < 2e-16 *** ## speciesChinstrap 9.908762 0.355289 27.889 < 2e-16 *** ## speciesGentoo 3.539179 0.499814 7.081 8.71e-12 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.419 on 329 degrees of freedom ## Multiple R-squared: 0.806, Adjusted R-squared: 0.8043 ## F-statistic: 455.8 on 3 and 329 DF, p-value: < 2.2e-16
modele.avec.interaction <- lm(bill_length_mm ~ body_mass_g * species, data = penguins) summary(modele.avec.interaction) #equation pour le manchot Gentoo?
## ## Call: ## lm(formula = bill_length_mm ~ body_mass_g * species, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.4164 -1.6630 0.0683 1.5444 9.3138 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 27.1129056 1.6324193 16.609 < 2e-16 *** ## body_mass_g 0.0031599 0.0004371 7.228 3.48e-12 *** ## speciesChinstrap 5.0612873 3.3101846 1.529 0.127 ## speciesGentoo -0.5750259 2.7941191 -0.206 0.837 ## body_mass_g:speciesChinstrap 0.0013028 0.0008832 1.475 0.141 ## body_mass_g:speciesGentoo 0.0009698 0.0006225 1.558 0.120 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.414 on 327 degrees of freedom ## Multiple R-squared: 0.8081, Adjusted R-squared: 0.8051 ## F-statistic: 275.3 on 5 and 327 DF, p-value: < 2.2e-16
En ajoutant une troisième variable à notre modèle, on arrive à un modèle de régression multiple. Attention: les variables indépendantes ne doivent pas être trop corrélées les unes aux autres (problème de la multicolinéarité).
modele <- lm(bill_length_mm ~ body_mass_g + species + year, data = penguins) summary(modele)
## ## Call: ## lm(formula = bill_length_mm ~ body_mass_g + species + year, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.4997 -1.6727 0.0599 1.4895 9.5913 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -5.978e+02 3.270e+02 -1.828 0.0685 . ## body_mass_g 3.751e-03 2.879e-04 13.031 < 2e-16 *** ## speciesChinstrap 9.935e+00 3.541e-01 28.053 < 2e-16 *** ## speciesGentoo 3.540e+00 4.978e-01 7.110 7.28e-12 *** ## year 3.101e-01 1.629e-01 1.904 0.0578 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.41 on 328 degrees of freedom ## Multiple R-squared: 0.8082, Adjusted R-squared: 0.8058 ## F-statistic: 345.5 on 4 and 328 DF, p-value: < 2.2e-16
On valide le modèle de la même façon.
par(mfrow=c(1,2)) plot(modele, which = 1:2)
La représentation graphique d’un modèle de régression multiple est plus difficile. Souvent, nous nous contentons de montrer les coefficients du modèle sur un forest plot.
sjPlot::plot_model(modele)
La représentation graphique d’un modèle de régression multiple est plus difficile. Souvent, nous nous contentons de montrer les coefficients du modèle sur un forest plot. Pour pouvoir comparer la taille des effets, il faut centrer-réduire toutes les variables quantitatives:
modele <- lm(bill_length_mm ~ scale(body_mass_g) + species + scale(year), data = penguins) sjPlot::plot_model(modele)
Quelques ressources utiles: